/*
 * =====================================================================================
 *
 *       Filename:  vm1.cpp
 *
 *    Description: Thomas algoritm
 *
 *        Version:  1.0
 *       Revision:  none
 *       Compiler:  gcc
 *
 *         Author:  Sergey Gvozdarev,
 *   Organization:  MIPT 016
 *
 *  =====================================================================================
 *
 */
#include <stdio.h>
#include <malloc.h>
#include <math.h>

int N;

double * x;

double * a;
double * b;
double * c;

/*
 * || -b0  c0    0 ...          ||
 * ||  a1 -b1   c1   0 ...      ||
 * ||   0  a2  -b2  c2   0 ...  ||
 * ||                           ||
 * ||                           ||
 * ||                           ||
 *
 */

void solve( ) 
{
    c[ 0] = c[ 0] / b[ 0];
    x[ 0] = x[ 0] / b[ 0];

    for ( int in = 1; in < N; in++) 
    {
        double m = 1.0 / ( b[ in] - a[ in] * c[ in - 1]);
        c[ in] = c[ in] * m;
        x[ in] = ( x[ in] - a[ in] * x[ in - 1]) * m;
    }

    for ( int in = N - 1; in > 0; in--)//?
    {
        x[ in] = x[ in] - c[ in] * x[ in + 1];
    }
}

double g( double x)
{
    return 1;
}

double h( double x)
{
    return x;
}

double s( double x)
{
    return -1;
}
double f( double x)
{
    return ( x * x + 1)/x;
}


int main( int argc, char** argv)
{
    
    // init

    double X_start = 0.5;
    double X_end   = 1.0;
    double tau     = 0.01;
    

    N = (int)( ( X_end - X_start) / tau);

    printf( "N: %d\n", N);
    a = (double*)malloc( sizeof( double) * ( N+1));
    b = (double*)malloc( sizeof( double) * ( N+1));
    c = (double*)malloc( sizeof( double) * ( N+1));
    x = (double*)malloc( sizeof( double) * ( N+1));

    for ( int i = 0; i < N + 1; i++)
    {
        a[ i] = b[ i] = c[ i] = x[ i] = 0.0;
    }


    tau = ( X_end - X_start) / (double)N;
    double tau_2 = tau * tau;
    double d1_tau   = 1 / tau;
    double d1_tau_2 = 1 / tau_2;

    b[ 0] = 1.0;
    x[ 0] = 0.0;

    b[ N-1] =    d1_tau;
    a[ N-1] =  - d1_tau;
    x[ N-1] = log( 2.0) + 1.0;

    double X_current = X_start;
    for ( int i = 1; i < N-1; i++)
    {
        a[ i] =    g( X_current - tau) * d1_tau_2 - 0.5 * h( X_current) * d1_tau;
        c[ i] =    g( X_current + tau) * d1_tau_2 + 0.5 * h( X_current) * d1_tau;
        b[ i] = - ( a[ i] + c[ i] - s( X_current));
        
        x[ i] = f( X_current);

        X_current = X_start + tau * i;
    }

    solve( );

    for ( int i = 0; i < N; i++)
    {
        printf( "%f\n", x[ i]);
    }

    free( a);
    free( b);
    free( c);
    free( x);
}
